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Excitonic couplings between (bacterio)chlorophyll molecules are necessary for simulating energy 
transport in photosynthetic complexes. Many techniques for calculating the couplings are in use, 
from the simple (but inaccurate) point-dipole approximation to fully quantum-chemical methods. 
We compared several approximations to determine their range of applicability, noting that the 
propagation of experimental uncertainties poses a fundamental limit on the achievable accuracy. 
In particular, the uncertainty in crystallographic coordinates yields an uncertainty of about 20% 
in the calculated couplings. Because quantum-chemical corrections are smaller than 20% in most 
biologically relevant cases, their considerable computational cost is rarely justihed. We therefore 
recommend the electrostatic TrEsp method across the entire range of molecular separations and 
orientations because its cost is minimal and it generally agrees with quantum-chemical calculations to 
better than the geometric uncertainty. We also caution against computationally optimizing a crystal 
structure before calculating couplings, as it can lead to large, uncontrollable errors. Understanding 
the unavoidable uncertainties can guard against striving for unrealistic precision; at the same time, 
detailed benchmarks can allow important qualitative questions—which do not depend on the precise 
values of the simulation parameters—to be addressed with greater confidence about the conclusions. 


Photosynthesis begins when (bacterio)chlorophyll 
molecules in an antenna complexes absorb light [1], 
creating molecular excited states. This triggers exci¬ 
tonic energy transfer (EET) [2], the migration of the 
excited-state energy through the network of (bacte- 
rio)chlorophyll until it decays or reaches the photosyn¬ 
thetic reaction center. 

The most common approach to modelling EET is to 
assume that each (bacterio) chlorophyll is a site—which 
can be in either the ground or excited states—so that 
the transfer of the excitation from one site to another is 
mediated by the coupling between them. In principle, 
the coupling can be calculated if the electronic structures 
and the relative positions of the two molecules are known. 
In practice, a full quantum-chemical treatment is often 
too expensive, which has led to various approximate 
methods for calculating excitonic couplings. This work is 
about determining the accuracy of these approximations 
and their range of applicability when they are applied 
to aggregates of bacteriochlorophylls. 

At separations that are much larger than the molecular 
dimensions, the leading term in the excitonic coupling is 
the dipole-dipole interaction of the two transition dipoles. 
This motivates the point-dipole (PD) approximation, 
which truncates all the higher-order contributions. Be¬ 
cause of its simplicity, the PD approximation has been 
widely used, even for calculating nearest-neighbor cou¬ 
plings where the small distance between the molecules 
would indicate that the approximation is inappropriate. 

Various methods—discussed below—go beyond the 
PD approximation and include more information about 
the electronic structure, but in a way that keeps the 
computation tractable. By comparing their performance 
with accurate quantum-chemical calculations, one might 
expect that a hierarchy of methods with defined distance 
cutoffs could be developed, so that the suitable method 
would be known for any particular intermolecular dis¬ 
tance. The difficulty with that approach is that setting 
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cutoffs depends on a subjective opinion about what is a 
tolerable error. 

However, an objective error threshold follows from the 
propagation of errors in the experimentally measured 
values that enter EET calculations. It is sufficient to 
use less accurate methods if the corrections would be 
drowned out by unavoidable uncertainties already in 
play. A particular inescapable error in EET calculations 
arises from the uncertainties in the atomic positions 
within each pigment. We show that a substantial uncer¬ 
tainty in the coupling is obtained by propagating the 
geometric uncertainties even if high-resolution crystal 
structures are used. In most cases, this uncertainty 
exceeds the quantum-mechanical corrections to the cou¬ 
pling, meaning that the substantially cheaper classical 
calculations are equally reliable. 


I. THEORY OF EXCITONIC COUPLINGS 

As two molecules are brought together, the Coulomb 
interaction of electrons and nuclei in one molecule with 
those in the other increases, meaning that the eigenstates 
of the isolated molecules are not eigenstates of the full 
Hamiltonian. However, if the intermolecular interaction 
is weak, it is often useful to think of the system as two 
interacting molecules as opposed to one supermolecule. 
To do so, one expands the full Hamiltonian—all the 
interactions between particles in either molecule—in the 
basis of molecular states, and the off-diagonal elements 
of that expansion are the couplings [2]. 

Excitonic couplings are interactions between excited 
states localized on different molecules. In the single- 
exciton manifold relevant to weak illumination, the 
Frenkel Hamiltonian of a system of two interacting two- 
level molecules is 


H = 




( 1 ) 


where J is the coupling between the donor and the 
acceptor, whose excitation energies are Ep, and Ea- 
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Method 


Short-range 

coupling 

Representation of 
transition density 

Cost 

Benefit 

PD 

Point dipole 

No 

Point dipole 

< 1 ms“ 

Lowest cost 

ED 

Extended dipole 

No 

Two charges 

< 1 ms“ 



TrEsp 

Transition charges from 
electrostatic potentials 

No 

Charges on atoms 

1 ms“ 



TDC 

Transition density cubes 

No 

Charges on a grid 

hours 



FED 

Fragment excitation difference 

Yes 

Full 

20 h 

Highest accuracy 


TABLE 1. Hierarchy of methods for calculating excitonic couplings. The cost is the typical time required to compute the 
coupling between two BChl molecules on a single processor. “ One-off tasks—the electronic-structure calculation of the 
transition density and, for TrEsp, fitting the atomic charges—are not included. ** Based on previous work [3]. 


The coupling is often described as containing short- 
range and long-range contributions. Short-range ef¬ 
fects include exchange, overlap of donor and acceptor 
wavefunctions, and exciton transfer mediated by charge- 
transfer states [4-7]. Because they depend on the spa¬ 
tial overlap between donor and acceptor wavefunctions, 
short-range terms decrease exponentially with distance 
and are consequently often neglected. 

Neglecting the short-range couplings leaves only the 
long-range Coulomb interaction, 


dcoui = 


1 


dTreo^r 


dro dr A 


P?g{^D)ptg{^A) 

\rD-y^A\ 


( 2 ) 


where and va are spatial coordinates and = 

if^{v)ip^{r) is the transition density between the ground 
and excited states of molecule X (either 

the donor D or the acceptor A). We return to the 
appropriate choice of the relative permittivity Sr below. 

Because Eq. 2 resembles the interaction of two charge 
distributions, the integral can be simplified using elemen¬ 
tary approximations from electrostatics. The various 
approximations—summarized in Table 1—differ in how 
they condense all the information in the continuous tran¬ 
sition densities into something more manageable and 
discrete. 

The simplest approximation, useful for intermolecular 
separations much larger than the sizes of the molecules, 
is the point-dipole approximation (PD), obtained as the 
lowest-order term in the multipole expansion of Eq. 2 [2], 


•/pD 


1 / dod4 _ ydc.r)(d,.r) -' 

dTrepEr \ J 


where r is the separation between the molecules 
(i.e., the centers of their transition densities) and 
djg and are their transition dipole moments, 
dx = JdrxTx pfgirx). 

However, PD is not accurate even at intermediate 
molecular separations [8-14]. A sequence of more ac¬ 
curate approximations can be obtained by representing 
the transition density as originating from an array of 
suitably chosen transition charges, 


</Coul 


1 

ATTSoEr 


EE 

iED jEA 





(4) 


The charges qi and their positions are chosen once 
and for all by htting them to the ab initio transition 


density, allowing each subsequent coupling calculation 
to be much faster. The different methods that have been 
used—extended dipole, TrEsp, and TDC—differ only in 
the number of transition charges used. 

In the extended dipole approximation (ED), two tran¬ 
sition charges are used per molecule, so that the Coulomb 
interaction becomes 



where are the positions of the positive and negative 
charges and ±5 is the magnitude of the charges. To be 
consistent with the point-dipole approximation at large 
separations, the charge 6 and the distance rj^ — must 
be chosen so that (5(rJ^ — r^^^) = djf, essentially making 
S an additional free parameter, whose tunability ensures 
ED agrees with the exact results better than PD. 

The opposite extreme is the transition density 
cube (TDC) method [15], where the transition charges 
are located on a Cartesian grid, making the method a 
direct numerical integration of Eq. 2. The main difficulty 
is that the grid needs to be Hne [3, 11, 16], containing 
many charges even in areas where the transition den¬ 
sity is negligible. For bacteriochlorophylls, the method 
converges when the number of charges in each molecule 
is around 500 000 [3], making the calculation of their 
pairwise interactions slow and prone to rounding errors. 

Between the two extremes of ED and TDC lies 
the transition monopole approximation (TMA) [17-19], 
which assigns one transition charge to each atom. If 
hydrogens are excluded, this leads to about 50 charges 
per BChl, a happy medium between 2 and 500 000. The 
atomic transition charges were initially assigned using a 
Mulliken or Hirshfeld population analysis [8], but this 
approach is not uniquely defined and was found to not ac¬ 
curately reproduce the shape of the transition density [3]. 
The method of transition charges from electrostatic po¬ 
tentials (TrEsp) [3, 20] solves this problem by htting 
the charges to best represent the transition density. For 
molecules with few atoms, additional htting parameters 
can be supplied by placing multipoles at each atom, 
which slightly improves the accuracy at the shortest 
separations [21, 22]. However, 50 charges offer enough 
free parameters that TrEsp is as accurate as TDC for 
chlorophylls [3], which is why we do not include TDC 
results in this study. 
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An alternative to transition charges and Eq. 4 is to 
expand the transition densities in a convenient chem¬ 
ical basis set and compute the resulting integrals in 
Eq. 2 using optimized quadrature techniques of quan¬ 
tum chemistry [23-27]. Like TDC, this approach gives 
the exact Coulomb coupling within the chosen basis, but 
is more expensive than TrEsp because it still requires 
numerical integration every time. As with TDC, we 
do not consider this approach here because TrEsp is 
sufficiently accurate. 

Once long-range couplings have been calculated using 
one of the methods surveyed above, any further im¬ 
provement must come from including short-range effects. 
Short-range couplings are particularly important in poly¬ 
cyclic aromatic hydrocarbons (PAH) [26, 28] and other 
flat molecules [8, 14, 29] whose lack of steric hindrance 
allows for tight packing. Among photosynthetic com¬ 
plexes, large short-range couplings have been reported in 
LH2 [30] and in the special pair of reaction centers [31]. 
We return to these cases below. 

The simplest situation is the coupling between a 
donor and an acceptor that are identical molecules. 
In that case, the eigenenergies of Eq. 1 are 
Ei ^2 = \{Ed + Ea) ± d, meaning that J can be ob¬ 
tained by halving the difference between the energies of 
the two excitonic states, 

J=i(Ai-A2), (6) 

which can be obtained from an electronic-structure cal¬ 
culation of the entire dimer. 

The same approach can be extended to het¬ 
erodimers [26, 30-33]. Scholes et al. used the eigen¬ 
states El 2 of H (from the quantum-chemical treatment 
of the dimer) and the site energies Eo and Ea (from 
the quantum-chemical treatment of the monomers) to 
calculate J [30]. Doing so assumed that the effective 
shifts in the site energies of the two molecules (induced 
by the presence of the other) are equal, an approxima¬ 
tion that was removed when the method was refined by 
Madjet et al [31]. They used the fact that heterodimer 
eigenstates are not fully delocalized, and that site-energy 
shifts can be computed from the extent of delocalization, 
which they obtained by comparing the monomer and 
dimer transition dipole moments. In this work, we use 
the closely related fragment excitation difference (FED) 
method [26, 33], which measures the delocalization in the 
eigenstates more directly, by determining the difference 
between the excitation densities on the two molecular 
fragments. Along with Ei^ 2 , this suffices to calculate J 
exactly, given a particular electronic-structure method 
and basis set. 


II. RESULTS AND DISCUSSION 

All coupling and energy calculations were performed 
on bacteriochlorophyll a (BChl a). In each case, the 
starting point was BChl a taken from the la position in 
the crystal structure of the LH2 complex of the purple 
bacterium Rhodopseudomonas acidophila [34]. 


A. Geometry optimization? 

An important preliminary question is which molecular 
geometry to use. It has been argued that the crystal 
structure is not reflective of the molecular configuration 
in vivo and that, consequently, the geometry should be 
computationally optimized before calculating the exci¬ 
tonic couplings [3, 31, 35, 36]. This optimization has 
generally been carried out using either HE or DFT. To 
the naked eye, the differences between the crystal struc¬ 
ture and the optimized geometry appear small (Figure 1) 
and one might think that this would result in only a 
minor correction to the couplings. 

However, the electronic properties of chlorophyll 
molecules are known to vary substantially with the 
molecular geometry [37-39]. To test the influence of 
geometry optimization on BChl couplings, we compared 
the electronic-structure results obtained with different 
geometry optimizations and with experiment. We cal¬ 
culated the Qj, transition energy at the crystal struc¬ 
ture geometry and following three different geometry 
optimizations: HF/6-31G* and DFT/B3LYP/6-31G* 
(gradient convergence criterion of 3 • 10“"^ Ah/oo in each 
case, Q-Chem 4.0 [40]), as well as molecular mechanics 
with Allinger’s MM3 forcefield [41] (gradient conver¬ 
gence criterion of 10“^ kcal/mol/A, TINKER 7.1 [42]). 
The transition energy itself was calculated using both 
CIS/6-31G* and TDDFT/B3LYP/6-31G*, as shown in 
Table 2. 

We found that the CIS calculation performed on 
the crystal structure gave the best agreement with 
the experimentally measured value for this transition, 
770 nm [43]. All geometry optimizations introduced er¬ 
rors much larger than the shifts in the Qy peak if the 
spectrum is measured in different solvents or following 
aggregation. This suggests that geometry optimization 
of this molecule prior to electronic-structure calcula¬ 
tions can introduce large and unnecessary errors. Con¬ 
sequently, all the following calculations were performed 
at the crystal-structure geometry using CIS/6-31G*. 

Table 2 also compares the Qy energies of the full 
BChl molecule and of the molecule with its phytyl tail 
removed. The removal of the tail changes the transition 
by at around 1 nm and can therefore be safely performed 
in the interest of computational speed. 


B. Coupling calculations 

Consistency is important to ensure a fair comparison 
of different methods. In particular, all the underly¬ 
ing electronic-structure calculations must use the same 
method and basis set, because different methods can give 
substantially different results [8, 12-14,16, 31, 32, 38, 44- 
46]. We used CIS/6-31G*, but this choice is not essential 
to our argument, since there is no reason to expect that 
other methods would be less affected by the geometric 
errors that are central to this paper. 

Two parameters in coupling calculations are difficult 
to determine ab initio: the relative permittivity Sr and 
the magnitude d of the transition dipole moment. The 
vacuum value = 1 is an appropriate choice for nearest- 
neighbor couplings in tightly packed aggregates, but 
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Optimised with 
HF/6-31G* 


FIG. 1. The crystal structure and the Hartree-Fock- 
optimized structure of BChl. Although the difference seems 
minor to the naked eye, it has a large influence on the tran¬ 
sition energies (see Table 2). 




Qy transition ( 

nm) 


With tail 

Tail removed 

Geometry 

CIS 

TDDFT 

CIS 

TDDFT 

Crystal structure 
Optimized: 

794 

621 

793 

621 

HF/6-31G* 

341 

559 

341 

559 

B3LYP/6-31G* 

664 

568 

667 

568 

MM3 

611 

578 

610 

578 

Experiment [43] 


770 



TABLE 2. The influence of geometry on the calculated 
wavelength of the Qy transition of BChl a. The transition was 
calculated using both CIS and TDDFT based on the crystal 
structure and on the geometries obtained by optimizing the 
crystal structure using different routines. The best agreement 
with experiment is obtained at the crystal structure geometry, 
indicating that geometry optimization may induce large and 
uncontrolled errors. In addition, removing the phytyl tail 
has a negligible effect, as expected. 


for more distant molecules, the Coulomb interaction is 
screened by the intervening medium. The appropriate 
choice of Sr for light-harvesting environments has been 
discussed extensively [35], and values from 1 to 2 have 
been used. This debate is beyond the scope of this work; 
instead, we report vacuum couplings, which can easily be 
adjusted by multiplication with l/e^. Our conclusions 
about the sizes of relative errors are unaffected, because 
all the couplings are scaled in the same proportion. 

The second uncertain parameter is the magnitude 
of the Qy transition dipole moment, which tends to 
be overestimated by CIS calculations [3, 15, 20, 30]. 
For the crystal-structure geometry, C1S/6-31G* predicts 
a Qy dipole moment of dciS = 10.45 D, significantly 
larger than the best experimental estimate of 6.1 D (in 
vacuum) [47]. Here again we make the simplest choice, 
reporting all results with the theoretically predicted 
value dciS i which can be corrected by multiplication with 
[d/dcisY, where d is the desired magnitude of the dipole 
moment [3, 15, 20, 30]. Again, the conclusions about 
the relative errors are unaffected because all couplings 
would be scaled by the same factor. 

Having made these preliminary choices, we compared 
the couplings predicted by PD, ED, TrEsp, and FED, 
calculated between two identical BChl molecules, dis¬ 
placed perpendicular to their bacteriochlorin rings by a 
separation ranging from 5 A to 20 A (Figure 2a). 

For PD, the magnitude and direction of the transition 
dipole moment were obtained from the CIS calculation. 
For ED, the two transition charges were chosen so that 
the ED coupling equaled the FED coupling at 20 A, 
giving a dipole length of 10.2 A. The dipole length is sen¬ 
sitive to the geometries used for the fitting, which is why 
our value differs from Renger’s 8.8 A [20]. For TrEsp, 
the best accuracy would be obtained by recalculating 
the transition charges at every geometry, but that would 
require an electronic-structure calculation and electro¬ 
static fitting each time, defeating TrEsp’s purpose as a 
fast method. Here, we use the original transition charges 
that the authors of TrEsp recommended be used even if 
the molecule undergoes slight configurational change [3]. 
However, because those charges were calculated for a pla¬ 
nar BChl a and predict a transition dipole of 10.11 D, we 


scaled all the charges by dcis/lO.ll D to ensure a consis¬ 
tent comparison with the other methods. For FED, we 
employed the routine implemented in Q-Chem 4.0 [40]. 

Figure 2b shows that all the coupling methods con¬ 
verge to a common value at large separations, as ex¬ 
pected. At the smallest separation in Figure 2b, the FED 
couplings differ from the TrEsp couplings by only 3%, 
indicating that the short-range contribution, ignored by 
TrEsp, is small compared to Jcoui- 

For separations under 5 A, FED results became diffi¬ 
cult to interpret because the molecules are so strongly 
coupled that it becomes impossible to speak of two 
coupled molecules and one must treat the dimer as a 
supermolecule. In particular, if the coupling becomes 
comparable to the spacing between electronic excited 
states, higher-energy transitions will contaminate the 
calculation and the two-state model of coupled Qy tran¬ 
sitions will fail. The breakdown of the two-state approx¬ 
imation is easily diagnosed in the parallel homodimer, 
because the approximation predicts that one of the two 
dimer states will be perfectly bright (twice the oscillator 
strength of the monomer) and the other perfectly dark 
(zero oscillator strength). For parallel BChls, this con¬ 
dition fails to hold around 4 A, making it dangerous to 
calculate couplings by simply halving the energetic gap. 

An Mg-Mg separation of 5 A is small and can only 
occur in configurations close to parallel because of the 
size of the BChl molecules. For comparison, the most 
strongly coupled naturally occurring BChls (with a 
known crystal structure) are the special pairs of reaction 
centers, with a Mg-Mg distance of about 8 A [48]. FED 
can fail in some of these cases as well, as we discuss 
below. 


C. Uncertainties 

Most theoretical calculations include experimentally 
measured quantities at some point and are therefore sub¬ 
ject to the propagation of errors. In particular, crystal- 
structure atomic coordinates carry uncertainties due to 
thermal motion and the limited resolution of the in- 
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FIG. 2. Excitonic couplings for the Qy transition of in vac¬ 
uum. (a) The geometry was taken from the crystal structure 
of LH2 of Rps. acidophila [34], with the second molecule 
a copy of the first, displaced perpendicular to the bacteri- 
ochlorin ring, (b) Comparison of four different methods 
(log-log scale). The error bars on the TrEsp values indicate 
the uncertainty in the coupling (20% on average), propa¬ 
gated from the uncertainty in the crystal-structure atom 
coordinates. Since the uncertainty is much larger than the 
difference between TrEsp and the fully quantum FED, the 
additional computational cost of FED is not justified. Inset: 
The deviation of TrEsp and ED from FED at short distances. 


couplings. We found that, for the parallel geometry, the 
resulting uncertainty in the TrEsp couplings (one stan¬ 
dard deviation) was about 20% across the whole range 
of separations. Since this is larger than the difference 
between the TrEsp and the FED values, we conclude 
that the considerable computational cost of FED is not 
justified at any separation. 

Our results are consistent with those of Arago and 
Troisi, who used molecular dynamics simulations of an¬ 
thracene crystals to find that the thermal nuclear mo¬ 
tion caused large fluctuations in excitonic couplings [28]. 
Although their results were dominated by short-range 
couplings—which are substantial in PAHs—our results 
confirm the importance of geometric fluctuations even 
for long-range couplings. 

The errors bars in Figure 2b are only the lower bound 
on the uncertainty in the couplings both because real¬ 
istic B factors are probably larger and because we did 
not consider other sources of uncertainty. We used the 
B factors from the LH2 crystal structure taken at 100 K; 
at physiological temperature, the thermal motion of 
the atoms would be greater. In addition, some simula¬ 
tions indicate that B factors computed in the course of 
ordinary structure refinement may underestimate the 
magnitude of thermal fluctuations [50]. The other kinds 
of uncertainty that may lead to substantially larger error 
bars include uncertainties in the choice of the transition 
dipole moment or the relative permittivity (especially 
at small-to-intermediate separations where there is little 
medium between the molecules [51]). Because these ad¬ 
ditional uncertainties could only increase the error bars, 
including them would strengthen our argument that the 
agreement between TrEsp and FED is better than the 
unavoidable error in the calculations. 


D. Other orientations 


strument. Those uncertainties should be propagated 
through the entire calculation. 

Uncertainties in crystallographic coordinates are spec¬ 
ified using Debye-Waller factors, also known as tem¬ 
perature factors or B factors [49]. They are given for 
each atom in the standard PDB crystallographic format, 
and are proportional to the mean-squared fluctuations 
of atomic positions. 



(7) 


Hence, the standard uncertaint y of eac h Cartesian co¬ 
ordinate is Ui = \J (u2)/3 = \/bJ%tt^. For the BChl 
molecule from the crystal structure of LH2 [34], the av- 

erage B factor is 16 A , corresponding to an uncertainty 
of 0.45 A in each Cartesian coordinate. 

We propagated these geometric uncertainties through 
our TrEsp calculations, and the results are shown as 
error bars in Figure 2b. For each calculation, the atomic 
positions in both molecules were chosen randomly from 
a normal distribution centered at the atom’s reported 
coordinates and with the standard deviation determined 
from the B factor. This was repeated with a sample of 
1000 random geometries to give the distribution of TrEsp 


To ensure that the good agreement between TrEsp 
and FED in Figure 2b was not peculiar to the parallel 
arrangement of the two BChls, we repeated the calcu¬ 
lation at different relative positions and orientations 
of the two molecules. We chose five molecular separa¬ 
tions (Mg-Mg distance) and calculated the ED, TrEsp, 
and FED couplings for 50 random orientations at each 
separation (see Figure 3b). Relative orientations were 
excluded if the two molecules collided, i.e., if any atom 
from the first molecule was closer than 1.5 A from any 
atom in the second. 

To determine the accuracy of the two approximate 
methods, the error of the ED and TrEsp calculations, 
with respect to the FED couplings, was calculated at 
each orientation. The typical relative error at each 
molecular separation was defined as 


AJ, 


ED/TrEsp 


|</fED — "/eD/TtEi 


iSp I 


(•/fed) 


( 8 ) 


where the averaging (•) was carried out over all the 
relative orientations at that molecular separation. This 
error is shown in Figure 3b, together with the standard 
deviation of the errors, | Jfed - •/ed/TtEspI/ (./fed)- 
The typical error of the ED approximation was large— 
about 50% of the FED coupling—indicating that the 
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b AJxrEsp I I 

AJed i i 

Average geometric uncertainty- 



FIG. 3. Accuracy of the TrEsp and ED methods, (a) The 
separation r between the magnesium atoms was fixed and 
the position (spherical coordinates r, Q, (j>) and orientation 
(Euler angles a, J, 7 ) of the second molecule were randomized 
(cases with overlapping molecules were rejected), (b) Error 
of TrEsp and ED, with respect to FED, averaged over 50 
random orientations at each separation. The average error 
is shown, together with error bars at one standard deviation. 
The dashed black line is the typical error arising from the 
uncertainty in the molecular geometry (the error bars in 
Figure 2b). ED can lead to large errors, indicating that 
the relatively good agreement in Figure 2b was coincidental. 
By contrast, TrEsp performs better than the propagated 
geometric uncertainty in almost all cases, confirming it as a 
robust method. 


method is not reliable in general and that the apparently 
good agreement in Figure 2b occurred because the ED 
dipole length was fitted to the FED data for the parallel 
arrangement. 


By contrast, TrEsp maintained the accuracy from 
Eigure 2b across the random orientational ensemble, 
with an average error of less than 6% at each separation. 
Eor all separations, the geometric uncertainty is more 
than one standard deviation higher than the average 
error of TrEsp. This indicates that it is unnecessary to 
use FED to calculate couplings between BChl molecules, 
as the limit of accuracy is not the choice between TrEsp 
and FED, but the uncertainty in the atomic positions. 


Pair 

r (A) 

Coupling (cm 


Relative error 

Jed 

t/TrEsp 

JpED 

AJed 

A JxrEsp 

LH2 







Jlalp 

9.0 

1099 

704 ± 326 

814 

-f35% 

-14% 

Jl/32a 

9.0 

214 

592 ± 238 

781 

-73% 

-24% 

RC 







JPlP2 

7.6 

35 

363 ± 371 

643“ 

-95% “ 

-44% “ 


TABLE 3. Comparison of ED, TrEsp and FED for some 
of the most strongly conpled bacteriochlorophylls in photo¬ 
synthetic complexes. The couplings were calculated, based 
on the crystal structures, for the nearest-neighbors in the 
B850 subunit of the LH2 complex of Rhodopseudomonas 
acidophila [34] and for the special pair in the reaction cen¬ 
ter of Rhodobacter sphaeroides [48]. As above, we assumed 
the molecules to be in vacuum and, for each molecule, we 
assumed transition dipoles predicted by CIS; therefore, the 
couplings should be scaled for comparison with previous 
work, r is the Mg-Mg separation and TrEsp values also 
include the geometric uncertainty. “ The FED calculation is 
unreliable due to the failure of the two-level approximation 
(see text). 


E. Most difficult cases 

The means and standard deviations in Figure 3b in¬ 
dicate that the difference between TrEsp and EED is 
less than the geometric uncertainty in the vast majority 
of cases. However, focusing on the mean and standard 
deviation risks ignoring the outliers, which may indi¬ 
cate additional failures of TrEsp. Of particular concern 
are configurations where the two BChls are parallel but 
offset from each other. The offset can give rise to a 
relatively large Mg-Mg distance even though portions 
of the two molecules might be close to each other. The 
offset-parallel arrangement arises in some natural com¬ 
plexes, giving rise to some of the most strongly coupled 
BChls in nature. To check the applicability of TrEsp to 
those cases, we calculated the couplings between three 
such pairs, as shown in Table 3. 

For the two nearest-neighbor couplings in the LH2 
complex of purple bacteria, the error of ED with respect 
to EED is over 50% on average, confirming its poor 
performance seen in Figure 3b. The error of TrEsp is, 
as expected, larger than the average in Eigure 3b, at 
14% and 24% for Jiaip and Jip 2 a, respectively, a result 
consistent with previous work finding short-range correc¬ 
tions of 17% and 24% in the two cases [30]. However, the 
geometric errors are also larger, because small displace¬ 
ments of particular atoms can have an outsized influence 
on the coupling when those atoms are close together. 
In the two cases from LH2, the average geometric un¬ 
certainty in the TrEsp couplings is 43%, substantially 
more than the 19% error with respect to FED. Thus, we 
can again conclude that the uncertainty in the atomic 
coordinates is a larger source of error than excluding 
the short-range contributions. 

Special pairs in reaction centers are even more strongly 
coupled, and the results in Table 3 indicate that all the 
methods fail. In particular, EED shows a considerable 
admixture of higher excited states, indicating that it 
is inappropriate to model this case as two coupled Qj, 
transitions (Eq. 1). Madjet et al. [31] checked the consis- 


































7 


tency of the same calculation by comparing the values 
of the Ya and parameters of their theory; these values 
were unequal by a large margin, especially for certain 
electronic-structure methods, indicating the failure of 
the effective two-state Hamiltonian. The difference be¬ 
tween our coupling and that of Madjet et al. is probably 
caused by their optimization of the geometry. The fail¬ 
ure of FED indicates that the special pair cannot be 
considered as two coupled molecules but should be seen 
as one unit. This can also be diagnosed from the se¬ 
vere failure of TrEsp, whose geometric uncertainty is 
over 100%. 

Therefore, the geometric uncertainty indicates not 
only the inherent uncertainty of TrEsp, but also the 
applicability of the two-state approximation. In difficult 
cases such as those just discussed, the close separation be¬ 
tween certain atoms will increase the uncertainty above 
the typical 20%, indicating that the calculation should 
be viewed with suspicion. 

III. CONCLUSIONS 

We have presented a detailed comparison of meth¬ 
ods for excitonic coupling calculations of bacteriochloro- 
phylls, with a focus on the uncertanties that arise due 
to the uncertainties in the atomic positions. For TrEsp, 
these geometric uncertainties are much larger than the 
disagreement between TrEsp and the more accurate 
FED, indicating that the short-range contributions to 


the couplings are almost always a minor correction in 
comparison with the error bars, making the computa¬ 
tional cost of obtaining them (6-7 orders of magnitude 
more than for TrEsp) unjustified. Furthermore, TrEsp 
is clearly preferrable over other the PD and ED approx¬ 
imations, both of which lead to errors much larger than 
the geometric uncertainty. Therefore, we recommend 
the use of TrEsp for the calculation of excitonic cou¬ 
plings between bacteriochlorophylls at all separations 
and orientations, with the warning that particularly 
large geometric uncertainties may indicate a failure of 
the two-state approximation. 

Our error analysis can be extended to site energies 
and other components of EET simulations in order to 
determine the overall sensitivity to uncertainties in the 
experimental data. This will allow us to identify the 
theoretical predictions that do not depend sensitively 
on microscopic details and are thus more likely to apply 
to a wide range of pigment-protein complexes. 
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